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摘 要 :扩展 有 限 元 法 通过 在 间断 区 域 引入 富 集 函数 ,在 处 理 强 弱 不 连续 问题 上 较 有 限 元 法 有 极 大 
的 优势 。 本 研究 给 出 了 基于 扩展 有 限 元 法 的 双 材 料 界面 裂纹 位 移 逼 近 方 程 及 相互 作用 积分 的 数值 

离散 方法 和 单元 刚度 矩阵 的 积分 策略 ,材料 界面 弱 不 连续 性 用 改进 扩展 有 限 元 模拟 ,裂纹 贯穿 部 分 

用 强 不 连续 的 Heaviside 函数 模拟 ,裂纹 尖端 分 别 用 2 种 不 同 渐 近 裂纹 尖端 富 集 函 数 模 拟 , 用 Matlab 

编制 相应 的 扩展 有 限 元 程序 。 算 例 表 明 ,数值 模拟 名 onn 5 Re ag EC 
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Numerical simulation of bimaterial interface cracks using 
extended finite element method 


SU Yi,CHEN Qinyuan 


(School of Aeronautical Engineering, Zhengzhou University of Aeronautics ,450046 Zhengzhou , China ) 


Abstract; The extended finite element method ( XFEM) has a great advantage over the finite element 
method ( FEM) in dealing with strong and weak discontinuities by introducing enrichment function into 
discontinuities. In this paper, bimaterial interface craks have been simulated using extended finite element 
method ( XFEM ). The displacement approximation equation of bimaterial interface crack is given. The 
integral method of the interaction integral and the integral strategy of element stiffness matrix are 
presented. Material discontinuity has been modeled by modified extended finite element method , and the 
crack penetration is simulated by Heaviside function ,the crack tip is simulated by two different asymptotic 
crack tip enrichment functions, and the corresponding propagation extended finite element program is 
compiled by Matlab. Good agreement between the numerical results and the reference solutions for 
bimaterial interface cracks problems is realized. 
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随 着 机 械 、 航 空 航天 和 生物 医学 应 用 对 多 性 能 
需求 的 (磨损 、 腐 蚀 、 耐 热 性 和 韧性 ) 日 益 增 加 , 双 材 
料 发 展 已 走 在 前 列 。 双 材料 的 整体 力学 行为 取决 于 
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方法 把 裂纹 分 成 裂纹 尖端 段 和 裂纹 主体 段 , 其 中 用 
Heaviside 函数 描述 裂纹 主体 段 ,用 二 维 渐进 裂纹 尖 
端 位 移 场 描 述 裂 尖端 段 。 采 用 文献 19] 提出 的 材料 


界面 的 力学 性 能 和 断裂 /疲劳 行为 。 但 因 其 上 下 两 
部 分 材料 的 机 械 性 能 差异 ,使 得 材料 界面 比较 薄弱 。 
因此 在 实际 使 用 过 程 中 , 双 材 料 的 断裂 破坏 ,往往 表 
现在 界面 的 破坏 。 如 层 合板 和 薄膜 基底 界面 的 脱 秋 
等 都 是 双 材 料 界面 问题 。 应 力 强度 因子 是 表征 界面 
裂纹 是 否 发 展 及 协助 判断 扩展 方向 的 重要 断裂 参 
数 。 因 此 精确 的 计算 应 力 强 度 因子 尤为 重要 。 求 解 
双 材 料 界面 裂纹 应 力 强 度 因子 常用 的 方法 有 加 料 有 
限 元 法 "边界 元 法 .扩展 有 限 元 法 趾 等 数值 计 
算 方法 。1999 年 , BELYTSCHKO 团队 ”基于 单位 
分 解法 在 有 限 元 的 基础 上 提出 了 扩展 有 限 元 法 。 它 
通过 在 位 移 场 上 对 局 部 强 弱 不 连续 区 域 增加 富 集 函 
数 来 描述 ,使 得 所 使 用 的 网 格 独立 于 结构 内 部 物理 
界面 ,从 而 避免 了 裂纹 尖端 等 高 应 力 区 网 格 划分 的 
困难 ,在 模拟 裂纹 生长 过 程 中 也 无 需 对 网 格 重新 划 
分 。 扩 展 有 限 元 法 在 处 理 间断 、 夹 杂 等 强 弱 不 连续 
问题 上 的 优点 使 其 自 被 提出 来 以 后 ,在 国际 上 得 到 
了 快速 的 发 展 和 广泛 的 应 用 。SUKUMAR 等 用 扩 
展 有 限 元 法 在 模拟 孔洞 夹杂 问题 时 ,用 符号 距离 
富 集 函数 来 描述 材料 界面 问题 ,但 数值 计算 收敛 性 
不 好 。 文 献 [4] 用 XFEM 模拟 双 材 料 界面 裂纹 ,用 
跨越 界面 裂纹 尖端 的 渐 近 位 移 场 来 作为 裂纹 尖端 的 
12 个 宣 集 函数 ,但 裂纹 边界 仍然 需要 沿 着 材料 界 
面 。 本 人 研究 用 文献 [9] 提 出 的 界面 富 集 函数 来 模拟 
界面 裂纹 的 材料 界面 ,用 Heaviside 函数 作为 贯穿 裂 
纹 富 集 函 数 ,用 改进 的 裂纹 尖端 渐 近 位 移 场 基 函 数 
作为 裂 尖 富 集 函 数 ,减少 自由 度 , 提 高 计算 效率 , i] 
时 能 够 不 降低 计算 精度 。 为 了 更 好 地 选择 裂 尖 富 集 
函数 ,同时 采用 2 种 不 同 的 裂 尖 富 集 函 数 模拟 裂纹 
尖端 位 移 场 的 奇异 性 ,并 编制 XFEM 程序 对 界面 裂 
纹 进行 研究 ,为 了 验证 程序 的 有 效 性 ,将 结果 与 数值 
结果 进行 比较 ,从 而 验证 程序 的 精确 性 。 


1 XFEM 的 基本 原理 


1.1 位 移 到 近 方程 的 建立 

扩展 有 限 元 法 区 别 于 有 限 元 法 在 于 它 基 于 单元 
分 解 通过 在 位 移 场 上 增加 富 集 也 数 来 描述 不 连续 问 
题 。 本 人 研究 为 了 更 好 地 模拟 裂纹 ,采用 文献 [5 ] 的 


界面 富 集 函 数 摘 述 双 材 料 界面 。 双 材料 界面 裂纹 如 
图 1 所 示 , 扩 展 有 限 元 法 的 位 移 逼 近 方 程 可 以 表 
示 为 ， 

u'(x) = Y N(Gx)u, + M, N(x)H()a,; + 


XO) Debi Y NGOw.G)e, C) 


xp IN QOO) ,N(x) ,Ni(x) ,N(x) 为 常规 有 限 元 法 
JÉ ERA SH. (x) e (x) ,W(x) 分 别 为 裂纹 面 ,裂纹 尖 
端 , 双 材 料 界面 的 富 集 孙 数 ;K 为 所 有 节点 集合 ,Ky、 
Ki Kj 分 别 表示 为 裂纹 面 、 裂 纹 尖 端 、 材料 界面 的 增 
强 节点 集合 ;aj、bi、c, 则 是 由 于 增强 节点 而 增加 的 
自由 度 。 
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图 1 双 材 料 界面 裂纹 
Fig.1 Bimaterial interface crack 
AR EX RBS E S PROS 
H(x)  H(£,(x)) -HCE Qx;)) (2) 
> Ẹ>0 


+1 


MOS i "ilo (3) 
裂纹 尖端 的 富 集 函 数 为 
Pi(X) 2o(x) -p(X) (4) 
Q,7 | /rsin 2 reos rsingsin 2. rsinfcos £] 
(5) 


式 中 ,r.9 是 裂纹 尖端 局 部 坐标 系 中 的 极 坐 标 , 裂纹 
尖端 为 极点 。 
在 确保 精度 的 情况 下 ,减少 自由 度 ,提高 计算 效 
率 ,根据 文献 [10] 将 裂纹 尖端 富 集 函数 变 为 
9, = [Vrsin (0/2 + m/4)Yrsingsin( 0/2 + m/4)] 
(6) 
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根据 文献 L111] 将 裂纹 尖端 宣 集 孔 数 表示 为 


HI FFR 
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; 为 载 集 列 阵 。 单 元 层次 上 的 和 分 别 可 表示 为 


$,(X) 122, 7 Wreoselnr) e “sin(0/2+m/4), 
Vreos( elnr) esin (0/2 + 0/4) , 
Vreos( elnr) e*^singsin( 9/2 + m/4) , 
Vrsin( elnr) e “sin(0/2 € q/4) , 
Vrsin( elnr) e^sin( 9/2 X 0/4) , 
Vrsin( elnr) e*^sin8sin( 8/2 + 5/4) | 
(7) 
Yn 为 双 材 料 界面 函数 ,根据 文献 [7] 用 扩展 有 
限 元 法 首次 模拟 富 集 函 数 , 双 材 料 界面 的 富 集 函 


lx) = | YN GOL, (8) 


AP, én 为 结 点 m 在 双 材 料 界面 水 平 集 函 数 中 的 
fü. 但 是 由 于 混合 单元 (单元 中 有 的 结 点 富 集 加 强 ， 
有 的 结 点 没有 加 强 ) 的 存在 影响 了 收敛 性 。 

为 了 元 服 这 个 问题 ,文献 19] 提 出 了 如 下 的 富 
R KZ, BI 


Pn = EN) [£6 | 7| NNQODE, 
] (9) 

£x) =, min [x -x || siga fæ) - f) 
(10) 


XB ux) 20 为 裂纹 曲线 方程 。 
模型 的 加 强 方式 如 图 2 所 示 。 


图 2 界面 裂纹 的 加 强 结 点 


Fig.2 Enrichment nodes of interface cracks 


1.2 XFEM 离散 方程 的 建立 


将 式 (1) 代 入 虚 功 方程 ,导出 XFEM 的 控制 方 
程 为 

Kô =F (11) 

式 中 :6 X625 xx R PUE: K 为 结构 总 刚度 矩阵 ;下 


uu ua ub uc 
ky? ky ky ky 


y 


au aa ab ac 
ky ky ky kj 
k; 二 bu ba bb bc ( 1 2 ) 
k? kj k; k; 
cu ca cb cc 
ky ky ky k; 


k; -[(BD'DBjdQ (r,s = u,a,b,c) (13) 
单元 荷载 列 阵 ,为 ， 


FELS SE SESS (14) 
单元 荷载 力 矢量 分 量 表示 为 


f= [Nitar + | Nba (15) 

f= Í N,HtdT + LNHbdn (16) 
T, 

k= |. Near + | Nigba (17) 

fi = [Natar + [ Napban (18) 


xh, B; UB; B? B; 为 形 函数 的 导数 矩阵 ,分 别 为 
常规 应 变 矩 阵 


c pM. 300 NAM 
un -| E Q(q, qo 4. CP 
裂纹 贯穿 单元 附加 应 变 和 矩阵 
(NH) 。 0 (N). 1 
un =| 0 — (NH), 254 ey 


裂纹 尖端 单元 附加 应 变 和 矩阵 
(Nipi) a 0 (CNien) 。 0 
0 GN, v 0 (Nip, ).,, 
(Nipi), CNigi) 。 Nipi) y (Ni) x 
(21) 


[B;]= 


FAAA E ET R H KIA AE RE 
(Na), 0 (Na) i. 
LIE i Yh (22) 
0 (Nu), (Na), 


L3 ”裂纹 尖端 应 力 强 度 因子 计算 
相互 作用 积分 的 形式 为 


a 
(23) 


式 中 :5 是 科 罗 内 区 尔 (Kronecker) 记号 ;Ww 8; 
o> 分 别 为 辅助 位 移 .应变 和 应 力 。 
用 散 度 定理 将 式 (23 ) 的 线 积分 转换 为 面积 积 


分 , 即 


y 


aux aux aux 
OakEi Ój — GjU;, — Oi Ui Jar 
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nu: 

CN 
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aux aux aux 
O 4€ ik Ôj Z OjUj;| ~ Oi Ui,l q ,dA 


(24) 
式 中 ,g 在 裂纹 尖端 处 取 1 ,在 边沿 C 处 取 0, 如 图 3 
所 示 。 


图 3 积分 区 域 


Fig.3 Integral region 


应 力 强度 因子 和 相互 作用 积分 的 关系 可 以 表 
ANN 


Ptol Oe ed + 天 TKT ) (25) 
式 中 
E ， 平面 应 力 
-| 所 , 平面 应 变 (26) 
e 是 双 材料 的 振荡 指数 , 即 
EE 
e= 元 gl E (27) 


ZB, a AG 描 述 双 材料 间 的 弹性 不 匹配 , 且 

qup +1) -W (K; +1) 

Tu +1) +u (Ki +1) 

qu -1) -p(x -1) 

pa C +1) tua( Ki +1) 

式 中 :py 为 材料 ii=1,2) 的 剪 切 模 量 ;x; 为 Kolosov 
常数 。 


(28) 


(29) 


] +v, (30) 
3 -vb, 平 面 应 变 
式 (25 ) 中 分 别 令 KE =1,KY =0, 此 时 把 1 记 
X, RE KP =1,KY =0, 此 时 把 1 记 为 可 以 求 
得 真实 场 作用 下 的 Kl1 和 Ki。o; 和 ji 由 扩展 有 限 
元 算得 ,w” .se 关 和 六 由 辅助 位 移 场 求 出 。 


! - EA 
Js gon, LI 


i : t MZ SIZ 
SCOPI LEES , 下 半 和 平面 
(31) 


3 =v. 
«= ,平面 应 力 
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当 求解 K| 时 ,辅助 位 移 场 为 


= D * 268sinsin 
fi p (32) 
= — C -26sinbcosp 
当 求解 Ki 时 ,辅助 位 移 场 为 
= — C *28sin6cos 
fi p (33) 
,- — D +2ôsinfsing 
AF,8,o,C,D 定义 见 文献 [11] 
aux 1 aux aux i * 
€i TG tu) (4j 21,2) (34) 
由 辅助 位 移 函 数 可 求 得 
f r af, 
iso eq) 
r jf. Fo 
T -A(Br., a =A( Bf. + 


1 


; 4u,cosh( me)’ 
式 中 :4=( 


roinn EFN 
p cosh ( me) 
DLEE ADLI o THARE se“ 求 得 ， 


上 半 平面 


2 相互 作用 积分 的 数值 离散 


在 XFEM 中 为 了 很 好 的 运用 相互 作用 积分 , 式 
(24) 可 以 被 离散 为 


=J 


ezl p= 


a £540,)q,| ,ll „W, (36) 
式 中 :e 是 积分 区 域 4 的 单元 的 数量 ;p, 在 单元 中 积 
分 点 的 数量 ;1J1 为 雅克 比 行列 式 ;w, 为 积分 点 p 
的 权重 。 实 际 位 移 的 导数 为 
4 4 
uj, = 2. us 机 2 bul Nya » $, +N, È or) + 


> aiVi H 十 > Cs UN, rs * Nus) (37) 
jekr mekp 


对 于 双 材 料 界面 裂纹 模型 ,建立 刚度 矩阵 时 , 材 
料 属性 由 积分 点 处 的 材料 属性 确定 。 由 于 在 形 函 数 
的 变 分 和 矩阵 中 含有 不 连续 的 函数 ,处 理 被 裂纹 分 割 
单元 时 ,文献 19 ,12-13 ] 给 出 了 积分 策略 。 在 本 研究 
中 ,处理 常 规 单元 采用 二 阶 高 斯 积分 ,混合 单元 采用 
六 阶 高 斯 积分 ,裂纹 贯穿 单元 ,需要 细 分 成 子 三 角形 
单元 ,而 后 每 个 子 三 角形 单元 采用 三 阶 高 斯 积分 ,而 
裂纹 尖端 单元 , 细 分 成 子 三 角形 单元 ,每 个 子 三 角形 


Pe 


| ( a g "E hl 
gg y OH 
1 
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单元 采用 七 阶 高 斯 积分 ,材料 界面 单元 采用 采用 六 
阶 高 斯 积分 。 


3 数值 算 例 


3.1 双 材 料 板 中 心 界 面 裂 纹 


如 图 4 所 示 , 模 型 是 高 为 2L = 200 mm, 宽 为 
2W =100 mm 的 双 材 料 板 ,其 中 心 有 一 界面 裂纹 , 板 受 
到 拉 伸 荷载 o, =9. 8 MPa 的 作用 ,裂纹 长 为 a/W = 
(0.4;0.8) ,材料 泊 松 比 v, 2 v, =0.3 ,假设 双 材 料 板 
处 于 平面 应 力 状 态 , 材 料 1 的 弹性 模 量 E, = 
2.058 x 10 MPa,E,/E, =1 ~100。 用 将 计算 得 到 的 
SIFs 通过 除 以 K, =o Vma 进 行 无 量 纲 化 处 理 , 即 
F, = K/K, ,并 将 MIYAZAKIL 和 NACGASHIMAL5 结 
果 与 本 研究 的 计算 结果 进行 比较 。 


| HH | 


7 


图 4 双 材料 板 中 心 的 界面 裂纹 (平面 应 力 ) 


Fig.4 An interface crack in the center of 


a bimaterial plate( plane stress ) 


从 图 5 ~6 可 以 看 出 ,本 研究 用 式 (6) RIS CT) 
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1.9 
1.8 
I 
16 Uu. 
—ec alW—-0.4 NAGASHIMA 
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Fig.5 F, of a center interface crack 
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图 6 中 心 界面 裂纹 FM 


Fig.6 F, of a center interface crack 


3.2. 无 限 大 双 材 料 板 中 心 界面 裂纹 问题 


如 图 7 所 示 , 一 个 二 维 弹性 双 材 料 平板 上 边界 
和 下 边界 受 单 向 偏远 拉 伸 载荷 o =1 MPa 作用 ,左右 
边界 受到 水 平 位 移 约束 u 20 mm, 在 模型 中 间 有 长 
度 为 2a 的 中 心 有 裂纹 。 板 矿 才 为 2L x 2L, 2L = 
60mm;a = 1 mm; E, = 1 x 10° MPa; E,/E, = 
2~1000。 根 据 文献 [4] 的 研究, 当 L/a > 20 时 ,该 板 
被 认为 是 无 限 大 板 。v =v, =0.3 ,假设 板 处 于 平面 


作为 界面 裂纹 尖端 富 集 函数 计算 的 无 量 纲 应 力 强 度 
-F 5 MIYAZAKI!“ NAGASHIMAL5 结果 近似 。 
说 明 用 式 (6) 和 式 (7) 都 可 以 作为 双 材 料 界面 裂纹 
裂纹 尖端 富 集 函 数 , 相 比 之 下 , 式 (7) 结 果 更 优 ,这 
是 因为 式 (7) 能够 更 好 地 模拟 界面 裂纹 尖端 位 移 
场 。F 的 绝对 值 随 着 Ig E,/E,). 的 增 大 而 减 小 ,而 
的 绝对 值 随 之 增加 。F; 的 绝对 值 随 裂纹 长 度 的 增 
加 而 增加 。 


应 变 状态 。 计 算 裂 纹 尖 端 应 力 强度 因子 K;, 无 量 纲 
1t, F; -K/K,, Hif Ky =0 Vma, J4 SUKUMAR! 
和 RICE 结果 与 本 研究 的 计算 结果 进行 比较 。 

从 图 8 ~9 可 以 看 出 ,本 研究 用 式 (6) 和 式 (7) 
作为 界面 裂纹 尖端 富 集 函数 计算 的 无 量 纲 应 力 强 度 
因子 与 参考 文献 结果 近似 ,了 的 易 合 效果 一 般 较 
丰 会 好 点 ,同样 用 式 (7 ) 模 拟 的 结果 更 好 一 些 。 但 
是 在 精度 要 求 不 高 的 情况 下 ,可 以 选择 式 (6)。 随 
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着 Ig(E,/E,) 的 增 大 ,了 的 绝对 值 随 之 增加 。 
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图 7 双 材 料 板 内 界面 裂纹 (平面 应 变 ) 
Fig.7 An interface crack in the center of 


a bimaterial plate( plane strain ) 
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图 9 ”中心 界面 裂纹 F, WEE 


Fig.9 F, of a center interface crack 
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4 结 it 

本 研究 介绍 了 扩展 有 限 元 在 双 材料 界面 裂纹 上 

的 应 用 。 在 求解 界面 柳 纹 的 应 力 强度 因子 时 用 界面 


富 集 函数 描述 界面 , 橡 尖 位 移 场 用 两 种 不 同 富 集 函 
数 模拟 。 数 值 模拟 的 结果 与 参考 文献 符合 地 较 好 ， 
但 当选 用 6 个 裂纹 尖端 富 集 函 数 精度 更 加 高 ,这 是 
因为 选用 的 裂纹 尖端 富 集 函 数 是 由 双 材料 界面 裂纹 
裂纹 尖端 渐 近 位 移 场 的 解析 解 的 基 郴 数 在 不 损失 精 


度 的 条 件 下 变换 的 富 集 函 数 ,但 6 个 富 


函数 计算 


量 大 ,而 用 2 个 计算 量 小 ,精度 稍 有 降低 。 
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